{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W3D1_RealNeurons/student/W3D1_Tutorial2.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy: Week 2, Day 3, Tutorial 2\n",
    "# Real Neurons: Effects of Input Correlation\n",
    "__Content creators:__ Qinglong Gu, Songtin Li, John Murray, Richard Naud, Arvind Kumar\n",
    "\n",
    "__Content reviewers:__ Maryam Vaziri-Pashkam, Ella Batty, Lorenzo Fontolan, Richard Gao, Matthew Krause, Spiros Chavlis, Michael Waskom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Tutorial Objectives\n",
    "In this tutorial, we will use the leaky integrate-and-fire (LIF) neuron model (see Tutorial 1) to study how they transform input correlations to output properties (transfer of correlations). In particular, we are going to write a few lines of code to:\n",
    "\n",
    "- inject correlated GWN in a pair of neurons\n",
    "\n",
    "- measure correlations between the spiking activity of the two neurons\n",
    "\n",
    "- study how the transfer of correlation depends on the statistics of the input, i.e. mean and standard deviation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:48.596870Z",
     "iopub.status.busy": "2021-05-25T01:12:48.595916Z",
     "iopub.status.idle": "2021-05-25T01:12:48.902441Z",
     "shell.execute_reply": "2021-05-25T01:12:48.903232Z"
    }
   },
   "outputs": [],
   "source": [
    "# Import libraries\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:48.909998Z",
     "iopub.status.busy": "2021-05-25T01:12:48.908286Z",
     "iopub.status.idle": "2021-05-25T01:12:48.991473Z",
     "shell.execute_reply": "2021-05-25T01:12:48.990944Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Figure Settings\n",
    "import ipywidgets as widgets       # interactive display\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "# use NMA plot style\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n",
    "my_layout = widgets.Layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.011428Z",
     "iopub.status.busy": "2021-05-25T01:12:49.004628Z",
     "iopub.status.idle": "2021-05-25T01:12:49.017008Z",
     "shell.execute_reply": "2021-05-25T01:12:49.017432Z"
    }
   },
   "outputs": [],
   "source": [
    "# @title Helper functions\n",
    "def default_pars(**kwargs):\n",
    "  pars = {}\n",
    "\n",
    "  ### typical neuron parameters###\n",
    "  pars['V_th'] = -55.     # spike threshold [mV]\n",
    "  pars['V_reset'] = -75.  # reset potential [mV]\n",
    "  pars['tau_m'] = 10.     # membrane time constant [ms]\n",
    "  pars['g_L'] = 10.       # leak conductance [nS]\n",
    "  pars['V_init'] = -75.   # initial potential [mV]\n",
    "  pars['V_L'] = -75.      # leak reversal potential [mV]\n",
    "  pars['tref'] = 2.       # refractory time (ms)\n",
    "\n",
    "  ### simulation parameters ###\n",
    "  pars['T'] = 400. # Total duration of simulation [ms]\n",
    "  pars['dt'] = .1  # Simulation time step [ms]\n",
    "\n",
    "  ### external parameters if any ###\n",
    "  for k in kwargs:\n",
    "    pars[k] = kwargs[k]\n",
    "\n",
    "  pars['range_t'] = np.arange(0, pars['T'], pars['dt'])  # Vector of discretized\n",
    "                                                         # time points [ms]\n",
    "  return pars\n",
    "\n",
    "\n",
    "def run_LIF(pars, Iinj):\n",
    "  \"\"\"\n",
    "  Simulate the LIF dynamics with external input current\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    Iinj       : input current [pA]. The injected current here can be a value or an array\n",
    "\n",
    "  Returns:\n",
    "    rec_spikes : spike times\n",
    "    rec_v      : mebrane potential\n",
    "  \"\"\"\n",
    "\n",
    "  # Set parameters\n",
    "  V_th, V_reset = pars['V_th'], pars['V_reset']\n",
    "  tau_m, g_L = pars['tau_m'], pars['g_L']\n",
    "  V_init, V_L = pars['V_init'], pars['V_L']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tref = pars['tref']\n",
    "\n",
    "  # Initialize voltage and current\n",
    "  v = np.zeros(Lt)\n",
    "  v[0] = V_init\n",
    "  Iinj = Iinj * np.ones(Lt)\n",
    "  tr = 0.\n",
    "\n",
    "  # simulate the LIF dynamics\n",
    "  rec_spikes = []   # record spike times\n",
    "  for it in range(Lt - 1):\n",
    "    if tr > 0:\n",
    "      v[it] = V_reset\n",
    "      tr = tr - 1\n",
    "    elif v[it] >= V_th:  # reset voltage and record spike event\n",
    "      rec_spikes.append(it)\n",
    "      v[it] = V_reset\n",
    "      tr = tref / dt\n",
    "\n",
    "    # calculate the increment of the membrane potential\n",
    "    dv = (-(v[it] - V_L) + Iinj[it] / g_L) * (dt / tau_m)\n",
    "\n",
    "    # update the membrane potential\n",
    "    v[it + 1] = v[it] + dv\n",
    "\n",
    "  rec_spikes = np.array(rec_spikes) * dt\n",
    "\n",
    "  return v, rec_spikes\n",
    "\n",
    "\n",
    "def my_GWN(pars, sig, myseed=False):\n",
    "  \"\"\"\n",
    "  Function that calculates Gaussian white noise inputs\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    mu         : noise baseline (mean)\n",
    "    sig        : noise amplitute (standard deviation)\n",
    "    myseed     : random seed. int or boolean\n",
    "                 the same seed will give the same random number sequence\n",
    "\n",
    "  Returns:\n",
    "    I          : Gaussian white noise input\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # Set random seed. You can fix the seed of the random number generator so\n",
    "  # that the results are reliable however, when you want to generate multiple\n",
    "  # realization make sure that you change the seed for each new realization\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # generate GWN\n",
    "  # we divide here by 1000 to convert units to sec.\n",
    "  I_GWN = sig * np.random.randn(Lt) * np.sqrt(pars['tau_m'] / dt)\n",
    "\n",
    "  return I_GWN\n",
    "\n",
    "\n",
    "def Poisson_generator(pars, rate, n, myseed=False):\n",
    "  \"\"\"\n",
    "  Generates poisson trains\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    rate       : noise amplitute [Hz]\n",
    "    n          : number of Poisson trains\n",
    "    myseed     : random seed. int or boolean\n",
    "\n",
    "  Returns:\n",
    "    pre_spike_train : spike train matrix, ith row represents whether\n",
    "                      there is a spike in ith spike train over time\n",
    "                      (1 if spike, 0 otherwise)\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve simulation parameters\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # set random seed\n",
    "  if myseed:\n",
    "      np.random.seed(seed=myseed)\n",
    "  else:\n",
    "      np.random.seed()\n",
    "\n",
    "  # generate uniformly distributed random variables\n",
    "  u_rand = np.random.rand(n, Lt)\n",
    "\n",
    "  # generate Poisson train\n",
    "  poisson_train = 1. * (u_rand < rate * (dt / 1000.))\n",
    "\n",
    "  return poisson_train\n",
    "\n",
    "def example_plot_myCC():\n",
    "  pars = default_pars(T=50000, dt=.1)\n",
    "\n",
    "  c = np.arange(10) * 0.1\n",
    "  r12 = np.zeros(10)\n",
    "  for i in range(10):\n",
    "    I1gL, I2gL = correlate_input(pars, mu=20.0, sig=7.5, c=c[i])\n",
    "    r12[i] = my_CC(I1gL, I2gL)\n",
    "\n",
    "  plt.figure()\n",
    "  plt.plot(c, r12, 'bo', alpha=0.7, label='Simulation', zorder=2)\n",
    "  plt.plot([-0.05, 0.95], [-0.05, 0.95], 'k--', label='y=x',\n",
    "           dashes=(2, 2), zorder=1)\n",
    "  plt.xlabel('True CC')\n",
    "  plt.ylabel('Sample CC')\n",
    "  plt.legend(loc='best')\n",
    "\n",
    "def LIF_output_cc(pars, mu, sig, c, bin_size, n_trials=20):\n",
    "  \"\"\" Simulates two LIF neurons with correlated input and computes output correlation\n",
    "\n",
    "  Args:\n",
    "  pars       : parameter dictionary\n",
    "  mu         : noise baseline (mean)\n",
    "  sig        : noise amplitute (standard deviation)\n",
    "  c          : correlation coefficient ~[0, 1]\n",
    "  bin_size   : bin size used for time series\n",
    "  n_trials   : total simulation trials\n",
    "\n",
    "  Returns:\n",
    "  r          : output corr. coe.\n",
    "  sp_rate    : spike rate\n",
    "  sp1        : spike times of neuron 1 in the last trial\n",
    "  sp2        : spike times of neuron 2 in the last trial\n",
    "  \"\"\"\n",
    "\n",
    "  r12 = np.zeros(n_trials)\n",
    "  sp_rate = np.zeros(n_trials)\n",
    "  for i_trial in range(n_trials):\n",
    "    I1gL, I2gL = correlate_input(pars, mu, sig, c)\n",
    "    _, sp1 = run_LIF(pars, pars['g_L'] * I1gL)\n",
    "    _, sp2 = run_LIF(pars, pars['g_L'] * I2gL)\n",
    "\n",
    "    my_bin = np.arange(0, pars['T'], bin_size)\n",
    "\n",
    "    sp1_count, _ = np.histogram(sp1, bins=my_bin)\n",
    "    sp2_count, _ = np.histogram(sp2, bins=my_bin)\n",
    "\n",
    "    r12[i_trial] = my_CC(sp1_count[::20], sp2_count[::20])\n",
    "    sp_rate[i_trial] = len(sp1) / pars['T'] * 1000.\n",
    "\n",
    "  return r12.mean(), sp_rate.mean(), sp1, sp2\n",
    "\n",
    "\n",
    "def plot_c_r_LIF(c, r, mycolor, mylabel):\n",
    "  z = np.polyfit(c, r, deg=1)\n",
    "  c_range = np.array([c.min() - 0.05, c.max() + 0.05])\n",
    "  plt.plot(c, r, 'o', color=mycolor, alpha=0.7, label=mylabel, zorder=2)\n",
    "  plt.plot(c_range, z[0] * c_range + z[1], color=mycolor, zorder=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "The helper function contains the:\n",
    "\n",
    "- Parameter dictionary: `default_pars( **kwargs)`\n",
    "- LIF simulator: `run_LIF`\n",
    "- Gaussian white noise generator: `my_GWN(pars, sig, myseed=False)`\n",
    "- Poisson type spike train generator: `Poisson_generator(pars, rate, n, myseed=False)`\n",
    "- Two LIF neurons with correlated inputs simulator: `LIF_output_cc(pars, mu, sig, c, bin_size, n_trials=20)`\n",
    "- Some additional plotting utilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1: Correlations (Synchrony)\n",
    "Correlation or synchrony in neuronal activity can be described for any readout of brain activity. Here, we are concerned with the spiking activity of neurons. \n",
    "\n",
    "In the simplest way, correlation/synchrony refers to coincident spiking of neurons, i.e., when two neurons spike together, they are firing in **synchrony** or are **correlated**. Neurons can be synchronous in their instantaneous activity, i.e., they spike together with some probability. However, it is also possible that spiking of a neuron at time $t$ is correlated with the spikes of another neuron with a delay (time-delayed synchrony). \n",
    "\n",
    "## Origin of synchronous neuronal activity:\n",
    "- Common inputs, i.e., two neurons are receiving input from the same sources. The degree of correlation of the shared inputs is proportional to their output correlation.\n",
    "- Pooling from the same sources. Neurons do not share the same input neurons but are receiving inputs from neurons which themselves are correlated.\n",
    "- Neurons are connected to each other (uni- or bi-directionally): This will only give rise to time-delayed synchrony. Neurons could also be connected via gap-junctions.\n",
    "- Neurons have similar parameters and initial conditions.\n",
    "\n",
    "## Implications of synchrony\n",
    "When neurons spike together, they can have a stronger impact on downstream neurons. Synapses in the brain are sensitive to the temporal correlations (i.e., delay) between pre- and postsynaptic activity, and this, in turn, can lead to the formation of functional neuronal networks - the basis of unsupervised learning (we will study some of these concepts in a forthcoming tutorial).\n",
    "\n",
    "Synchrony implies a reduction in the dimensionality of the system. In addition, correlations, in many cases, can impair the decoding of neuronal activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.025510Z",
     "iopub.status.busy": "2021-05-25T01:12:49.024637Z",
     "iopub.status.idle": "2021-05-25T01:12:49.082168Z",
     "shell.execute_reply": "2021-05-25T01:12:49.081709Z"
    },
    "outputId": "be989be4-76f0-4d3d-9d18-51b3d4c5c5a7"
   },
   "outputs": [],
   "source": [
    "# @title Video 1: Input & output correlations\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"nsAYFBcAkes\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## How to study the emergence of correlations\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "A simple model to study the emergence of correlations is to inject common inputs to a pair of neurons and measure the output correlation as a function of the fraction of common inputs. \n",
    "\n",
    "Here, we are going to investigate the transfer of correlations by computing the correlation coefficient of spike trains recorded from two unconnected LIF neurons, which received correlated inputs.\n",
    "\n",
    "\n",
    "The input current to LIF neuron $i$ $(i=1,2)$ is:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{I_i}{g_L} =\\mu_i + \\sigma_i (\\sqrt{1-c}\\xi_i + \\sqrt{c}\\xi_c) \\quad (1)\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mu_i$ is the temporal average of the current. The Gaussian white noise $\\xi_i$ is independent for each neuron, while $\\xi_c$ is common to all neurons. The variable $c$ ($0\\le c\\le1$) controls the fraction of common and independent inputs. $\\sigma_i$ shows the variance of the total input.\n",
    "\n",
    "So, first, we will generate correlated inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 247
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.088340Z",
     "iopub.status.busy": "2021-05-25T01:12:49.087776Z",
     "iopub.status.idle": "2021-05-25T01:12:49.092641Z",
     "shell.execute_reply": "2021-05-25T01:12:49.091861Z"
    },
    "outputId": "eb01ad76-9757-4384-d008-f6fde57cafed"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "#@markdown Execute this cell to get a function for generating correlated GWN inputs\n",
    "def correlate_input(pars, mu=20., sig=7.5, c=0.3):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    mu         : noise baseline (mean)\n",
    "    sig        : noise amplitute (standard deviation)\n",
    "    c.         : correlation coefficient ~[0, 1]\n",
    "\n",
    "  Returns:\n",
    "    I1gL, I2gL : two correlated inputs with corr. coe. c\n",
    "  \"\"\"\n",
    "\n",
    "  # generate Gaussian whute noise xi_1, xi_2, xi_c\n",
    "  xi_1 = my_GWN(pars, sig)\n",
    "  xi_2 = my_GWN(pars, sig)\n",
    "  xi_c = my_GWN(pars, sig)\n",
    "\n",
    "  # Generate two correlated inputs by Equation. (1)\n",
    "  I1gL = mu + np.sqrt(1. - c) * xi_1 + np.sqrt(c) * xi_c\n",
    "  I2gL = mu + np.sqrt(1. - c) * xi_2 + np.sqrt(c) * xi_c\n",
    "\n",
    "  return I1gL, I2gL\n",
    "\n",
    "print(help(correlate_input))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 1: Compute the correlation\n",
    "\n",
    "The _sample correlation coefficient_ between two input currents $I_i$ and $I_j$ is defined as the sample covariance of $I_i$ and $I_j$ divided by the square root of the sample variance of $I_i$ multiplied with the square root of the sample variance of $I_j$. In equation form:  \n",
    "\n",
    "\\begin{align}\n",
    "r_{ij} &= \\frac{cov(I_i, I_j)}{\\sqrt{var(I_i)} \\sqrt{var(I_j)}}\\\\\n",
    "cov(I_i, I_j) &= \\sum_{k=1}^L (I_i^k -\\bar{I}_i)(I_j^k -\\bar{I}_j) \\\\\n",
    "var(I_i) &= \\sum_{k=1}^L (I_i^k -\\bar{I}_i)^2\n",
    "\\end{align}\n",
    "\n",
    "where $\\bar{I}_i$ is the sample mean, k is the time bin, and L is the length of $I$.  This means that $I_i^k$ is current i at time $k\\cdot dt$. Note that the equations above are not accurate for sample covariances and variances as they should be additionally divided by L-1 - we have dropped this term because it cancels out in the sample correlation coefficient formula.\n",
    "\n",
    "The _sample correlation coefficient_ may also be referred to as the _sample Pearson correlation coefficient_. Here, is a beautiful paper that explains multiple ways to calculate and understand correlations [Rodgers and Nicewander 1988](https://www.stat.berkeley.edu/~rabbee/correlation.pdf).\n",
    "\n",
    "In this exercise, we will create a function, `my_CC` to compute the sample correlation coefficient between two time series. Note that while we introduced this computation here in the context of input currents, the sample correlation coefficient is used to compute the correlation between any two time series - we will use it later on binned spike trains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.096957Z",
     "iopub.status.busy": "2021-05-25T01:12:49.096419Z",
     "iopub.status.idle": "2021-05-25T01:12:49.099976Z",
     "shell.execute_reply": "2021-05-25T01:12:49.099465Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_CC(i, j):\n",
    "  \"\"\"\n",
    "  Args:\n",
    "    i, j  : two time series with the same length\n",
    "\n",
    "  Returns:\n",
    "    rij   : correlation coefficient\n",
    "  \"\"\"\n",
    "  ########################################################################\n",
    "  ## TODO for students: compute rxy, then remove the NotImplementedError #\n",
    "  # Tip1: array([a1, a2, a3])*array([b1, b2, b3]) = array([a1*b1, a2*b2, a3*b3])\n",
    "  # Tip2: np.sum(array([a1, a2, a3])) = a1+a2+a3\n",
    "  # Tip3: square root, np.sqrt()\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student exercise: compute the sample correlation coefficient\")\n",
    "  ########################################################################\n",
    "  # Calculate the covariance of i and j\n",
    "  cov = ...\n",
    "  # Calculate the variance of i\n",
    "  var_i = ...\n",
    "  # Calculate the variance of j\n",
    "  var_j = ...\n",
    "  # Calculate the correlation coefficient\n",
    "  rij = ...\n",
    "\n",
    "  return rij\n",
    "\n",
    "\n",
    "# Uncomment the line after completing the my_CC function\n",
    "# example_plot_myCC()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 431
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.105818Z",
     "iopub.status.busy": "2021-05-25T01:12:49.105231Z",
     "iopub.status.idle": "2021-05-25T01:12:49.979629Z",
     "shell.execute_reply": "2021-05-25T01:12:49.978790Z"
    },
    "outputId": "b5fa5ff8-d328-46d4-c6aa-fa4dcd95bb7a"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D3_RealNeurons/solutions/W2D3_Tutorial2_Solution_03e44bdc.py)\n",
    "\n",
    "*Example output:*\n",
    "\n",
    "<img alt='Solution hint' align='left' width=558 height=413 src=https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/tutorials/W2D3_RealNeurons/static/W2D3_Tutorial2_Solution_03e44bdc_0.png>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 2: Measure the correlation between spike trains\n",
    "\n",
    "After recording the spike times of the two neurons, how can we estimate their correlation coefficient? \n",
    "\n",
    "In order to find this, we need to bin the spike times and obtain two time series. Each data point in the time series is the number of spikes in the corresponding time bin. You can use `np.histogram()` to bin the spike times.\n",
    "\n",
    "Complete the code below to bin the spike times and calculate the correlation coefficient for two Poisson spike trains. Note that `c` here is the ground-truth correlation coefficient that we define."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 247
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.986466Z",
     "iopub.status.busy": "2021-05-25T01:12:49.985873Z",
     "iopub.status.idle": "2021-05-25T01:12:49.988078Z",
     "shell.execute_reply": "2021-05-25T01:12:49.988558Z"
    },
    "outputId": "ae08f9d2-1388-4940-ac1f-40c7992497b0"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Execute this cell to get a function for generating correlated Poisson inputs (generate_corr_Poisson)\n",
    "\n",
    "\n",
    "def generate_corr_Poisson(pars, poi_rate, c, myseed=False):\n",
    "  \"\"\"\n",
    "  function to generate correlated Poisson type spike trains\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    poi_rate   : rate of the Poisson train\n",
    "    c.         : correlation coefficient ~[0, 1]\n",
    "\n",
    "  Returns:\n",
    "    sp1, sp2   : two correlated spike time trains with corr. coe. c\n",
    "  \"\"\"\n",
    "\n",
    "  range_t = pars['range_t']\n",
    "\n",
    "  mother_rate = poi_rate / c\n",
    "  mother_spike_train = Poisson_generator(pars, rate=mother_rate,\n",
    "                                         n=1, myseed=myseed)[0]\n",
    "  sp_mother = range_t[mother_spike_train > 0]\n",
    "\n",
    "  L_sp_mother = len(sp_mother)\n",
    "  sp_mother_id = np.arange(L_sp_mother)\n",
    "  L_sp_corr = int(L_sp_mother * c)\n",
    "\n",
    "  np.random.shuffle(sp_mother_id)\n",
    "  sp1 = np.sort(sp_mother[sp_mother_id[:L_sp_corr]])\n",
    "\n",
    "  np.random.shuffle(sp_mother_id)\n",
    "  sp2 = np.sort(sp_mother[sp_mother_id[:L_sp_corr]])\n",
    "\n",
    "  return sp1, sp2\n",
    "\n",
    "\n",
    "print(help(generate_corr_Poisson))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:49.996206Z",
     "iopub.status.busy": "2021-05-25T01:12:49.994881Z",
     "iopub.status.idle": "2021-05-25T01:12:49.996793Z",
     "shell.execute_reply": "2021-05-25T01:12:49.997211Z"
    }
   },
   "outputs": [],
   "source": [
    "def corr_coeff_pairs(pars, rate, c, trials, bins):\n",
    "  \"\"\"\n",
    "  Calculate the correlation coefficient of two spike trains, for different\n",
    "  realizations\n",
    "\n",
    "  Args:\n",
    "      pars   : parameter dictionary\n",
    "      rate   : rate of poisson inputs\n",
    "      c      : correlation coefficient ~ [0, 1]\n",
    "      trials  : number of realizations\n",
    "      bins   : vector with bins for time discretization\n",
    "\n",
    "  Returns:\n",
    "    r12      : correlation coefficient of a pair of inputs\n",
    "  \"\"\"\n",
    "\n",
    "  r12 = np.zeros(n_trials)\n",
    "\n",
    "  for i in range(n_trials):\n",
    "    ##############################################################\n",
    "    ## TODO for students: Use np.histogram to bin the spike time #\n",
    "    ##   e.g., sp1_count, _= np.histogram(...)\n",
    "    # Use my_CC() compute corr coe, compare with c\n",
    "    # Note that you can run multiple realizations and compute their r_12(diff_trials)\n",
    "    # with the defined function above. The average r_12 over trials can get close to c.\n",
    "    # Note: change seed to generate different input per trial\n",
    "    # Fill out function and remove\n",
    "    raise NotImplementedError(\"Student exercise: compute the correlation coefficient\")\n",
    "    ##############################################################\n",
    "\n",
    "    # Generate correlated Poisson inputs\n",
    "    sp1, sp2 = generate_corr_Poisson(pars, ..., ..., myseed=2020+i)\n",
    "\n",
    "    # Bin the spike times of the first input\n",
    "    sp1_count, _ = np.histogram(..., bins=...)\n",
    "\n",
    "    # Bin the spike times of the second input\n",
    "    sp2_count, _ = np.histogram(..., bins=...)\n",
    "\n",
    "    # Calculate the correlation coefficient\n",
    "    r12[i] = my_CC(..., ...)\n",
    "\n",
    "  return r12\n",
    "\n",
    "\n",
    "poi_rate = 20.\n",
    "c = 0.2  # set true correlation\n",
    "pars = default_pars(T=10000)\n",
    "\n",
    "# bin the spike time\n",
    "bin_size = 20  # [ms]\n",
    "my_bin = np.arange(0, pars['T'], bin_size)\n",
    "n_trials = 100  # 100 realizations\n",
    "\n",
    "# Uncomment to test your function\n",
    "# r12 = corr_coeff_pairs(pars, rate=poi_rate, c=c, trials=n_trials, bins=my_bin)\n",
    "# print(f'True corr coe = {c:.3f}')\n",
    "# print(f'Simu corr coe = {r12.mean():.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Sample output\n",
    "\n",
    "```\n",
    "True corr coe = 0.200\n",
    "Simu corr coe = 0.197\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 52
    },
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:50.004578Z",
     "iopub.status.busy": "2021-05-25T01:12:50.004011Z",
     "iopub.status.idle": "2021-05-25T01:12:50.151867Z",
     "shell.execute_reply": "2021-05-25T01:12:50.152305Z"
    },
    "outputId": "572dcd51-b621-4ab1-8ecd-a3f34ded05cc"
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D3_RealNeurons/solutions/W2D3_Tutorial2_Solution_e5eaac3e.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Investigate the effect of input correlation on the output correlation\n",
    "\n",
    "Now let's combine the aforementioned two procedures. We first generate the correlated inputs by Equation (1). Then we inject the correlated inputs $I_1, I_2$ into a pair of neurons and record their output spike times. We continue measuring the correlation between the output and \n",
    "investigate the relationship between the input correlation and the output correlation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Drive a neuron with correlated inputs and visualize its output\n",
    "In the following, you will inject correlated GWN in two neurons. You need to define the mean (`gwn_mean`), standard deviation (`gwn_std`), and input correlations (`c_in`).\n",
    "\n",
    "We will simulate $10$ trials to get a better estimate of the output correlation. Change the values in the following cell for the above variables (and then run the next cell) to explore how they impact the output correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:50.158036Z",
     "iopub.status.busy": "2021-05-25T01:12:50.157035Z",
     "iopub.status.idle": "2021-05-25T01:12:50.159007Z",
     "shell.execute_reply": "2021-05-25T01:12:50.159398Z"
    }
   },
   "outputs": [],
   "source": [
    "# Play around with these parameters\n",
    "\n",
    "pars = default_pars(T=80000, dt=1.)  # get the parameters\n",
    "c_in = 0.3  # set input correlation value\n",
    "gwn_mean = 10.\n",
    "gwn_std = 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 484
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:50.166361Z",
     "iopub.status.busy": "2021-05-25T01:12:50.165791Z",
     "iopub.status.idle": "2021-05-25T01:12:52.405752Z",
     "shell.execute_reply": "2021-05-25T01:12:52.405295Z"
    },
    "outputId": "c4b6b9cd-47fc-4505-b5de-ce5c6546aefc"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Do not forget to execute this cell to simulate the LIF\n",
    "\n",
    "\n",
    "bin_size = 10.  # ms\n",
    "\n",
    "starttime = time.perf_counter()  # time clock\n",
    "\n",
    "r12_ss, sp_ss, sp1, sp2 = LIF_output_cc(pars, mu=gwn_mean, sig=gwn_std, c=c_in,\n",
    "                                        bin_size=bin_size, n_trials=10)\n",
    "\n",
    "# just the time counter\n",
    "endtime = time.perf_counter()\n",
    "timecost = (endtime - starttime) / 60.\n",
    "print(f\"Simulation time = {timecost:.2f} min\")\n",
    "\n",
    "print(f\"Input correlation = {c_in}\")\n",
    "print(f\"Output correlation = {r12_ss}\")\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(sp1, np.ones(len(sp1)) * 1, '|', ms=20, label='neuron 1')\n",
    "plt.plot(sp2, np.ones(len(sp2)) * 1.1, '|', ms=20, label='neuron 2')\n",
    "plt.xlabel('time (ms)')\n",
    "plt.ylabel('neuron id.')\n",
    "plt.xlim(1000, 8000)\n",
    "plt.ylim(0.9, 1.2)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Think!\n",
    "- Is the output correlation always smaller than the input correlation? If yes, why?\n",
    "- Should there be a systematic relationship between input and output correlations? \n",
    "\n",
    "You will explore these questions in the next figure but try to develop your own intuitions first!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "Lets vary `c_in` and plot the relationship between the `c_in` and output correlation. This might take some time depending on the number of trials. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:12:52.413696Z",
     "iopub.status.busy": "2021-05-25T01:12:52.413127Z",
     "iopub.status.idle": "2021-05-25T01:13:12.639010Z",
     "shell.execute_reply": "2021-05-25T01:13:12.638049Z"
    },
    "outputId": "af94ba78-ed5a-4574-9344-b4fcdc101fc7"
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Don't forget to execute this cell!\n",
    "\n",
    "pars = default_pars(T=80000, dt=1.)  # get the parameters\n",
    "bin_size = 10.\n",
    "c_in = np.arange(0, 1.0, 0.1)  # set the range for input CC\n",
    "r12_ss = np.zeros(len(c_in))  # small mu, small sigma\n",
    "\n",
    "starttime = time.perf_counter() # time clock\n",
    "for ic in range(len(c_in)):\n",
    "  r12_ss[ic], sp_ss, sp1, sp2 = LIF_output_cc(pars, mu=10.0, sig=10.,\n",
    "                                              c=c_in[ic], bin_size=bin_size,\n",
    "                                              n_trials=10)\n",
    "\n",
    "endtime = time.perf_counter()\n",
    "timecost = (endtime - starttime) / 60.\n",
    "print(f\"Simulation time = {timecost:.2f} min\")\n",
    "\n",
    "plt.figure(figsize=(7, 6))\n",
    "plot_c_r_LIF(c_in, r12_ss, mycolor='b', mylabel='Output CC')\n",
    "plt.plot([c_in.min() - 0.05, c_in.max() + 0.05],\n",
    "         [c_in.min() - 0.05, c_in.max() + 0.05],\n",
    "         'k--', dashes=(2, 2), label='y=x')\n",
    "\n",
    "plt.xlabel('Input CC')\n",
    "plt.ylabel('Output CC')\n",
    "plt.legend(loc='best', fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:12.643552Z",
     "iopub.status.busy": "2021-05-25T01:13:12.642977Z",
     "iopub.status.idle": "2021-05-25T01:13:12.648210Z",
     "shell.execute_reply": "2021-05-25T01:13:12.647744Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D3_RealNeurons/solutions/W2D3_Tutorial2_Solution_71e76f4d.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 3: Correlation transfer function\n",
    "The above plot of input correlation vs. output correlation is called the __correlation transfer function__ of the neurons. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 3.1: How do the mean and standard deviation of the GWN affect the correlation transfer function?\n",
    "\n",
    "The correlations transfer function appears to be linear. The above can be taken as the input/output transfer function of LIF neurons for correlations, instead of the transfer function for input/output firing rates as we had discussed in the previous tutorial (i.e., F-I curve).\n",
    "\n",
    "What would you expect to happen to the slope of the correlation transfer function if you vary the mean and/or the standard deviation of the GWN?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:12.659717Z",
     "iopub.status.busy": "2021-05-25T01:13:12.658759Z",
     "iopub.status.idle": "2021-05-25T01:13:43.402457Z",
     "shell.execute_reply": "2021-05-25T01:13:43.401573Z"
    },
    "outputId": "6db28df6-3256-4517-ae67-1c187c3cdf2c"
   },
   "outputs": [],
   "source": [
    "#@markdown Execute this cell to visualize correlation transfer functions\n",
    "\n",
    "pars = default_pars(T=80000, dt=1.) # get the parameters\n",
    "no_trial = 10\n",
    "bin_size = 10.\n",
    "c_in = np.arange(0., 1., 0.2)  # set the range for input CC\n",
    "r12_ss = np.zeros(len(c_in))   # small mu, small sigma\n",
    "r12_ls = np.zeros(len(c_in))   # large mu, small sigma\n",
    "r12_sl = np.zeros(len(c_in))   # small mu, large sigma\n",
    "\n",
    "starttime = time.perf_counter()  # time clock\n",
    "for ic in range(len(c_in)):\n",
    "  r12_ss[ic], sp_ss, sp1, sp2 = LIF_output_cc(pars, mu=10.0, sig=10.,\n",
    "                                              c=c_in[ic], bin_size=bin_size,\n",
    "                                              n_trials=no_trial)\n",
    "  r12_ls[ic], sp_ls, sp1, sp2 = LIF_output_cc(pars, mu=18.0, sig=10.,\n",
    "                                              c=c_in[ic], bin_size=bin_size,\n",
    "                                              n_trials=no_trial)\n",
    "  r12_sl[ic], sp_sl, sp1, sp2 = LIF_output_cc(pars, mu=10.0, sig=20.,\n",
    "                                              c=c_in[ic], bin_size=bin_size,\n",
    "                                              n_trials=no_trial)\n",
    "endtime = time.perf_counter()\n",
    "timecost = (endtime - starttime) / 60.\n",
    "print(f\"Simulation time = {timecost:.2f} min\")\n",
    "\n",
    "plt.figure(figsize=(7, 6))\n",
    "plot_c_r_LIF(c_in, r12_ss, mycolor='b', mylabel=r'Small $\\mu$, small $\\sigma$')\n",
    "plot_c_r_LIF(c_in, r12_ls, mycolor='y', mylabel=r'Large $\\mu$, small $\\sigma$')\n",
    "plot_c_r_LIF(c_in, r12_sl, mycolor='r', mylabel=r'Small $\\mu$, large $\\sigma$')\n",
    "plt.plot([c_in.min() - 0.05, c_in.max() + 0.05],\n",
    "         [c_in.min() - 0.05, c_in.max() + 0.05],\n",
    "         'k--', dashes=(2, 2), label='y=x')\n",
    "plt.xlabel('Input CC')\n",
    "plt.ylabel('Output CC')\n",
    "plt.legend(loc='best', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Think!\n",
    "Why do both the mean and the standard deviation of the GWN affect the slope of the correlation transfer function? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:43.406665Z",
     "iopub.status.busy": "2021-05-25T01:13:43.406098Z",
     "iopub.status.idle": "2021-05-25T01:13:43.411305Z",
     "shell.execute_reply": "2021-05-25T01:13:43.410840Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D3_RealNeurons/solutions/W2D3_Tutorial2_Solution_2deb4ccb.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Section 3.2: What is the rationale behind varying $\\mu$ and $\\sigma$?\n",
    "The mean and the variance of the synaptic current depends on the spike rate of a Poisson process. We can use [Campbell's theorem](https://en.wikipedia.org/wiki/Campbell%27s_theorem_(probability)) to estimate the mean and the variance of the synaptic current:\n",
    "\n",
    "\\begin{align}\n",
    "\\mu_{\\rm syn} = \\lambda J \\int P(t) \\\\\n",
    "\\sigma_{\\rm syn} = \\lambda J \\int P(t)^2 dt\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $\\lambda$ is the firing rate of the Poisson input, $J$ the amplitude of the postsynaptic current and $P(t)$ is the shape of the postsynaptic current as a function of time. \n",
    "\n",
    "Therefore, when we varied $\\mu$ and/or $\\sigma$ of the GWN, we mimicked a change in the input firing rate. Note that, if we change the firing rate, both $\\mu$ and $\\sigma$ will change simultaneously, not independently. \n",
    "\n",
    "Here, since we observe an effect of $\\mu$ and $\\sigma$ on correlation transfer, this implies that the input rate has an impact on the correlation transfer function.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Think!\n",
    "\n",
    "- What are the factors that would make output correlations smaller than input correlations? (Notice that the colored lines are below the black dashed line)\n",
    "- What does it mean for the correlation in the network?\n",
    "- Here we have studied the transfer of correlations by injecting GWN. But in the previous tutorial, we mentioned that GWN is unphysiological. Indeed, neurons receive colored noise (i.e., Shot noise or OU process). How do these results obtained from injection of GWN apply to the case where correlated spiking inputs are injected in the two LIFs? Will the results be the same or different?\n",
    "\n",
    "Reference\n",
    "- De La Rocha, Jaime, et al. \"Correlation between neural spike trains increases with firing rate.\" Nature (2007) (https://www.nature.com/articles/nature06028/)\n",
    "\n",
    "- Bujan AF, Aertsen A, Kumar A. Role of input correlations in shaping the variability and noise correlations of evoked activity in the neocortex. Journal of Neuroscience. 2015 Jun 3;35(22):8611-25. (https://www.jneurosci.org/content/35/22/8611)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "text",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:43.416028Z",
     "iopub.status.busy": "2021-05-25T01:13:43.415466Z",
     "iopub.status.idle": "2021-05-25T01:13:43.418138Z",
     "shell.execute_reply": "2021-05-25T01:13:43.418575Z"
    }
   },
   "source": [
    "[*Click for solution*](https://github.com/NeuromatchAcademy/course-content/tree/master//tutorials/W2D3_RealNeurons/solutions/W2D3_Tutorial2_Solution_39d29f52.py)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "In this tutorial, we studied how the input correlation of two LIF neurons is mapped to their output correlation. Specifically, we:\n",
    "\n",
    "- injected correlated GWN in a pair of neurons,\n",
    "\n",
    "- measured correlations between the spiking activity of the two neurons, and\n",
    "\n",
    "- studied how the transfer of correlation depends on the statistics of the input, i.e., mean and standard deviation.\n",
    "\n",
    "Here, we were concerned with zero time lag correlation. For this reason, we restricted estimation of correlation to instantaneous correlations. If you are interested in time-lagged correlation, then we should estimate the cross-correlogram of the spike trains and find out the dominant peak and area under the peak to get an estimate of output correlations. \n",
    "\n",
    "We leave this as a future to-do for you if you are interested."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 1: Example of a conductance-based LIF model\n",
    "Above, we have written code to generate correlated Poisson spike trains. You can write code to stimulate the LIF neuron with such correlated spike trains and study the correlation transfer function for spiking input and compare it to the correlation transfer function obtained by injecting correlated GWNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 354
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:43.431718Z",
     "iopub.status.busy": "2021-05-25T01:13:43.431095Z",
     "iopub.status.idle": "2021-05-25T01:13:43.433708Z",
     "shell.execute_reply": "2021-05-25T01:13:43.434160Z"
    },
    "outputId": "7b16fb3f-864a-4ef5-a858-39b5b9d7d78e"
   },
   "outputs": [],
   "source": [
    "# @title Function to simulate conductance-based LIF\n",
    "\n",
    "def run_LIF_cond(pars, I_inj, pre_spike_train_ex, pre_spike_train_in):\n",
    "  \"\"\"\n",
    "  conductance-based LIF dynamics\n",
    "\n",
    "  Args:\n",
    "    pars               : parameter dictionary\n",
    "    I_inj              : injected current [pA]. The injected current here can\n",
    "                         be a value or an array\n",
    "    pre_spike_train_ex : spike train input from presynaptic excitatory neuron\n",
    "    pre_spike_train_in : spike train input from presynaptic inhibitory neuron\n",
    "\n",
    "  Returns:\n",
    "    rec_spikes         : spike times\n",
    "    rec_v              : mebrane potential\n",
    "    gE                 : postsynaptic excitatory conductance\n",
    "    gI                 : postsynaptic inhibitory conductance\n",
    "  \"\"\"\n",
    "\n",
    "  # Retrieve parameters\n",
    "  V_th, V_reset = pars['V_th'], pars['V_reset']\n",
    "  tau_m, g_L = pars['tau_m'], pars['g_L']\n",
    "  V_init, E_L = pars['V_init'], pars['E_L']\n",
    "  gE_bar, gI_bar = pars['gE_bar'], pars['gI_bar']\n",
    "  VE, VI = pars['VE'], pars['VI']\n",
    "  tau_syn_E, tau_syn_I = pars['tau_syn_E'], pars['tau_syn_I']\n",
    "  tref = pars['tref']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "\n",
    "  # Initialize\n",
    "  tr = 0.\n",
    "  v = np.zeros(Lt)\n",
    "  v[0] = V_init\n",
    "  gE = np.zeros(Lt)\n",
    "  gI = np.zeros(Lt)\n",
    "  Iinj = I_inj * np.ones(Lt)  # ensure I has length Lt\n",
    "\n",
    "  if pre_spike_train_ex.max() == 0:\n",
    "    pre_spike_train_ex_total = np.zeros(Lt)\n",
    "  else:\n",
    "    pre_spike_train_ex_total = pre_spike_train_ex * np.ones(Lt)\n",
    "\n",
    "  if pre_spike_train_in.max() == 0:\n",
    "    pre_spike_train_in_total = np.zeros(Lt)\n",
    "  else:\n",
    "    pre_spike_train_in_total = pre_spike_train_in * np.ones(Lt)\n",
    "\n",
    "  # simulation\n",
    "  rec_spikes = []  # recording spike times\n",
    "  for it in range(Lt - 1):\n",
    "    if tr > 0:\n",
    "      v[it] = V_reset\n",
    "      tr = tr - 1\n",
    "    elif v[it] >= V_th:   # reset voltage and record spike event\n",
    "      rec_spikes.append(it)\n",
    "      v[it] = V_reset\n",
    "      tr = tref / dt\n",
    "    # update the synaptic conductance\n",
    "    gE[it+1] = gE[it] - (dt / tau_syn_E) * gE[it] + gE_bar * pre_spike_train_ex_total[it + 1]\n",
    "    gI[it+1] = gI[it] - (dt / tau_syn_I) * gI[it] + gI_bar * pre_spike_train_in_total[it + 1]\n",
    "\n",
    "    # calculate the increment of the membrane potential\n",
    "    dv = (-(v[it] - E_L) - (gE[it + 1] / g_L) * (v[it] - VE) - \\\n",
    "            (gI[it + 1] / g_L) * (v[it] - VI) + Iinj[it] / g_L) * (dt / tau_m)\n",
    "\n",
    "    # update membrane potential\n",
    "    v[it + 1] = v[it] + dv\n",
    "\n",
    "  rec_spikes = np.array(rec_spikes) * dt\n",
    "\n",
    "  return v, rec_spikes, gE, gI\n",
    "\n",
    "\n",
    "print(help(run_LIF_cond))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Interactive Demo: Correlated spike input to an LIF neuron \n",
    "\n",
    "In the following you can explore what happens when the neurons receive correlated spiking input.\n",
    "\n",
    "You can vary the correlation between excitatory input spike trains. For simplicity, the correlation between inhibitory spike trains is set to 0.01.\n",
    "\n",
    "Vary both excitatory rate and correlation and see how the output correlation changes. Check if the results are qualitatively similar to what you observed previously when you varied the $\\mu$ and $\\sigma$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 634,
     "referenced_widgets": [
      "2a40d29dc8014f39ae26cd7a103ed121",
      "0e038564c631440cada9810cc24f3f45",
      "4239c3f494be4c298fb9777c3774ed27",
      "fe40759b66fa488b8d855c829298fa48",
      "bbcc5c90735b4c46b2c55721fee2dc15",
      "06a2d60056274ffbac0965da71e7cffe",
      "bd215fdd48314c409d45ccec15ca41c7",
      "b49eeb974089428b8c16f94ed8d6eed6",
      "db07d79071d94f43a2e15ad1819557a0",
      "0778a54dfe164c7ea3fdcadba3fb145b",
      "fc4e64907dbe4001a6a0568a1fee7d08"
     ]
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:43.480321Z",
     "iopub.status.busy": "2021-05-25T01:13:43.478585Z",
     "iopub.status.idle": "2021-05-25T01:13:44.191635Z",
     "shell.execute_reply": "2021-05-25T01:13:44.191109Z"
    },
    "outputId": "3ab2a5a2-769a-4388-c52d-7ef391702736"
   },
   "outputs": [],
   "source": [
    "# @title\n",
    "\n",
    "# @markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "my_layout.width = '450px'\n",
    "@widgets.interact(\n",
    "  pwc_ee=widgets.FloatSlider(0.3, min=0.05, max=0.99, step=0.01,\n",
    "                             layout=my_layout),\n",
    "  exc_rate=widgets.FloatSlider(1e3, min=500., max=5e3, step=50.,\n",
    "                               layout=my_layout),\n",
    "  inh_rate=widgets.FloatSlider(500., min=300., max=5e3, step=5.,\n",
    "                               layout=my_layout),\n",
    ")\n",
    "\n",
    "\n",
    "def EI_isi_regularity(pwc_ee, exc_rate, inh_rate):\n",
    "  pars = default_pars(T=1000.)\n",
    "  # Add parameters\n",
    "  pars['V_th'] = -55.     # spike threshold [mV]\n",
    "  pars['V_reset'] = -75.  # reset potential [mV]\n",
    "  pars['tau_m'] = 10.     # membrane time constant [ms]\n",
    "  pars['g_L'] = 10.       # leak conductance [nS]\n",
    "  pars['V_init'] = -65.   # initial potential [mV]\n",
    "  pars['E_L'] = -75.      # leak reversal potential [mV]\n",
    "  pars['tref'] = 2.       # refractory time (ms)\n",
    "\n",
    "  pars['gE_bar'] = 4.0    # [nS]\n",
    "  pars['VE'] = 0.         # [mV] excitatory reversal potential\n",
    "  pars['tau_syn_E'] = 2.  # [ms]\n",
    "  pars['gI_bar'] = 2.4    # [nS]\n",
    "  pars['VI'] = -80.       # [mV] inhibitory reversal potential\n",
    "  pars['tau_syn_I'] = 5.  # [ms]\n",
    "\n",
    "  my_bin = np.arange(0, pars['T']+pars['dt'], .1)  # 20 [ms] bin-size\n",
    "\n",
    "  # exc_rate = 1e3\n",
    "  # inh_rate = 0.4e3\n",
    "  # pwc_ee = 0.3\n",
    "  pwc_ii = 0.01\n",
    "  # generate two correlated spike trains for excitatory input\n",
    "  sp1e, sp2e = generate_corr_Poisson(pars, exc_rate, pwc_ee)\n",
    "\n",
    "  sp1_spike_train_ex, _ = np.histogram(sp1e, bins=my_bin)\n",
    "  sp2_spike_train_ex, _ = np.histogram(sp2e, bins=my_bin)\n",
    "\n",
    "  # generate two uncorrelated spike trains for inhibitory input\n",
    "  sp1i, sp2i = generate_corr_Poisson(pars, inh_rate, pwc_ii)\n",
    "  sp1_spike_train_in, _ = np.histogram(sp1i, bins=my_bin)\n",
    "  sp2_spike_train_in, _ = np.histogram(sp2i, bins=my_bin)\n",
    "\n",
    "  v1, rec_spikes1, gE, gI = run_LIF_cond(pars, 0, sp1_spike_train_ex, sp1_spike_train_in)\n",
    "  v2, rec_spikes2, gE, gI = run_LIF_cond(pars, 0, sp2_spike_train_ex, sp2_spike_train_in)\n",
    "\n",
    "  # bin the spike time\n",
    "  bin_size = 20  # [ms]\n",
    "  my_bin = np.arange(0, pars['T'], bin_size)\n",
    "  spk_1, _ = np.histogram(rec_spikes1, bins=my_bin)\n",
    "  spk_2, _ = np.histogram(rec_spikes2, bins=my_bin)\n",
    "\n",
    "  r12 = my_CC(spk_1, spk_2)\n",
    "  print(f\"Input correlation = {pwc_ee}\")\n",
    "  print(f\"Output correlation = {r12}\")\n",
    "\n",
    "  plt.figure(figsize=(14, 7))\n",
    "  plt.subplot(211)\n",
    "  plt.plot(sp1e, np.ones(len(sp1e)) * 1, '|', ms=20,\n",
    "           label='Exc. input 1')\n",
    "  plt.plot(sp2e, np.ones(len(sp2e)) * 1.1, '|', ms=20,\n",
    "           label='Exc. input 2')\n",
    "  plt.plot(sp1i, np.ones(len(sp1i)) * 1.3, '|k', ms=20,\n",
    "           label='Inh. input 1')\n",
    "  plt.plot(sp2i, np.ones(len(sp2i)) * 1.4, '|k', ms=20,\n",
    "           label='Inh. input 2')\n",
    "\n",
    "  plt.ylim(0.9, 1.5)\n",
    "  plt.legend()\n",
    "  plt.ylabel('neuron id.')\n",
    "\n",
    "  plt.subplot(212)\n",
    "  plt.plot(pars['range_t'], v1, label='neuron 1')\n",
    "  plt.plot(pars['range_t'], v2, label='neuron 2')\n",
    "  plt.xlabel('time (ms)')\n",
    "  plt.ylabel('membrane voltage $V_{m}$')\n",
    "  plt.tight_layout()\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "\n",
    "Above, we are estimating the output correlation for one trial. You can modify the code to get a trial average of output correlations.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Bonus 2: Ensemble Response\n",
    "\n",
    "Finally, there is a short BONUS lecture video on the firing response of an ensemble of neurons to time-varying input. There are no associated coding exercises - just enjoy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2021-05-25T01:13:44.204026Z",
     "iopub.status.busy": "2021-05-25T01:13:44.203460Z",
     "iopub.status.idle": "2021-05-25T01:13:44.238551Z",
     "shell.execute_reply": "2021-05-25T01:13:44.238978Z"
    },
    "outputId": "66a909fd-a0ca-42d2-b20c-b07acb74c2e9"
   },
   "outputs": [],
   "source": [
    "#@title Video 2 (Bonus): Response of ensemble of neurons to time-varying input\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"78_dWa4VOIo\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W2D3_Tutorial2",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
